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Abstract 

The indefinite sign of the Hamiltonian constraint means that so- 
lutions to Einstein's equations must achieve a delicate balance — often 
among numerically large terms that nearly cancel. If numerical errors 
cause a violation of the Hamiltonian constraint, the failure of the del- 
icate balance could lead to qualitatively wrong behavior rather than 
just decreased accuracy. This issue is different from instabilities caused 
by constraint-violating modes. Examples of stable numerical simula- 
tions of collapsing cosmological spacetimes exhibiting local mixmaster 
dynamics with and without Hamiltonian constraint enforcement are 
presented. 

1 Introduction 1 

It gives me great pleasure to dedicate this note to Mike Ryan in honor of 
his 60th birthday. It is especially rewarding to be able to include several of 
Mike's favorite research themes including Mixmaster dynamics and spatially 
inhomogeneous cosmologies — where he has made important contributions. 

The 3 + 1 formulation of Einstein's equations has served as the starting 
point for most numerical simulations. While the precise set of Einstein's 
equations most suitable for simulations is a subject of active investigation 
(for a review see [1]), the basic structure of the ADM form — evolution equa- 
tions for (something related to) the induced spatial metric and (something 

^This paper is based on a talk given at the 17th International Conference on General 
Relativity and Gravitation. 
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related to) the spatial extrinsic curvature plus (some form of) the Hamilto- 
nian and momentum constraints — will most likely remain. As is well known, 
the Einstein evolution equations preserve the constraints but the discretized 
evolution equations need not do so. To date, most simulations of binary 
black holes and/or neutron stars have solved only the evolution equations. 
Experience with these simulations has demonstrated that failure to solve the 
constraints allows the growth of constraint violating modes that can cause 
the codes to crash. The main argument against numerical enforcement of 
the constraints has been the cost in compute time of solving the elliptic 
constraint equations. Advances in computational power, improved elliptic 
solvers, and the growing awareness of the problem of instabilities caused by 
constraint violation have led to increased interest in constrained evolution 

[2]- 

In this note, I shall discuss a further reason to consider constrained 
evolution — especially for the Hamiltonian constraint. The Hamiltonian con- 
straint is given schematically by 

n° = -±= (k^k 13 - k 2 ) - Vf? R [h l3 } = o (i) 

where K i3 - and hy are respectively the extrinsic curvature and induced metric 
of the spatial hypersurface, h is the determinant of hij, 3 R[hij] is the scalar 
curvature of h^, and K is the trace of K^. Note first the indefinite sign 
of HP which can appear in either the "kinetic" terms involving K i3 - or in 
the "potential" term containing 3 R. Thus to maintain the solution 7i° = 
requires a delicate balance among the variables. Failure to achieve this 
balance in a numerical simulation can yield not only numerical instability 
but also — even in a stable evolution — qualitatively incorrect behavior. In 
the following, three examples of collapsing cosmological spacetimes with 3, 
2, and 1 spatial Killing field will be used to illustrate qualitatively incorrect 
behavior that results if the Hamiltonian constraint is allowed to evolve freely 
in a simulation. 

Collapsing cosmological spacetimes are characterized by two main types of 
behavior — asymptotic velocity term dominance (AVTD) and local Mixmaster 
dynamics (LMD) — first described by Belinskii et al (BKL) (e.g. [3]). Both 
types of behavior arise when the dynamics becomes local — variables at each 
spatial point evolve as a separate spatially homogeneous universe. The basic 
building block of the approach to the singularity is the Kasner universe, 
characterized by fixed collapse rates along the principle spatial axes. This 
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is in contrast to Bianchi IX (Mixmaster) collapse where an infinite sequence 
of "bounces" off the spatial scalar curvature change one Kasner epoch into 
another with different fixed collapse rates. The relationship between one 
Kasner epoch and the next includes bounces that are sensitive to initial 
conditions and can be determined from conservation of "momentum" through 
the bounce [3, 4]. 

Collapsing cosmologies exhibit LMD if, at (almost) 2 every spatial point, 
bounces from one Kasner-like epoch to another may be demonstrated and 
shown to obey the relevant bounce laws. In contrast, AVTD behavior is 
characterized by a final Kasner epoch at (almost) every spatial point after a 
possible final bounce. Rigorous demonstration of AVTD behavior has been 
provided [5]. With the exception of spatially homogeneous cosmologies [6], 
LMD has been plausibly demonstrated only through numerical simulation 
[7,8,9]. 

In the remainder of this note, we shall consider three examples where 
LMD can be seen only if the Hamiltonian constraint is explicitly enforced. 
These examples are the spatially homogeneous vacuum, diagonal Bianchi 
IX cosmology [3, 4, 10], the two spatial Killing field generic T 2 -symmetric 
vacuum spacetime [11, 8], and the one spatial Killing field vacuum U(l)- 
symmetric cosmology [12, 13]. While failure to enforce the Hamiltonian con- 
straint can lead to instability, we shall focus on simulations that evolve stably 
but indicate that the approach to the singularity is AVTD rather than LMD. 

2 Examples 

2.1 Bianchi IX (Mixmaster) Cosmology 

The Mixmaster universe is described by the metric [4] 

ds 2 = _ e 3n dr 2 + e 2n/ e 2/j\ d(J i d(T j (2) 

V / ij 

where = diag(— 2/9+, / 9 + +2v / 3 / 9_, , (3 + — 2-\/3/3-), ^, P± depend only on the 
time r, and the spatial 1-forms a 1 satisfy the appropriate SU (2) relationship 
for Bianchi IX. Einstein's equations may be found from the variation of the 

2 The word "almost" will be used in this context to indicate that set-of-measure-zero 
exceptions are known. Details of these can usually be found in the cited references. 



3 



superhamiltonian (lapse times the Hamiltonian constraint) specialized for 
these models as [4] 



2H° = = -pi + p\ + p 2 _ + v(n,p ± ) 



(3) 



where pn, p± are canonically conjugate to £l,(3± and 

V(Q,(3±) = e 4n_8/3 + + e 40 + 4 ' 3 ++ 4 v / 3/3- + e 4Q+4/3_+4x/3/3_ + 



(4) 



where the ellipsis indicates terms that are (almost) always negligible. In the 
absence of the potential V (with a 1 = dx l ), the Kasner solution is obtained. 
Eq. (3) defines the dynamics in minisuperspace (MSS). The Kasner solution 
represents the free particle in MSS. For the Kasner solution, Eq. (3) may be 
written as 



where v± = —p±/pn- The addition of the potential (4) causes (almost) every 
Kasner epoch to end in a bounce off one of the exponential terms in V. After 
the bounce, the behavior is again described by (5) but with different Kasner 
parameters. As first discussed by BKL, every Kasner epoch can be identiifed 
by a single parameter u (related to the anisotropic collapse rates) such that 
the n + 1 st Kasner epoch is related to the n th one in Mixmaster dynamics 
through 



As the Mixmaster universe evolves toward the singularity, the ratio of the 
duration of the bounce to the duration of the Kasner epoch goes to zero. The 
role of the Hamiltonian constraint (3) — as it becomes arbitrarily close to (5) 
except at the bounce — is to keep the configuration on the "Kasner circle" 
defined by K — 1. A typical evolution may be found in [14] , Fig. 5. 

A numerical simulation of Mixmaster collapse with a code that can follow 
the evolution through hundreds of bounces [10] shows qualitatively different 
behavior depending upon whether or not the Hamiltonian constraint is en- 
forced. Figure 1 shows a comparison of the Kasner circle indicator K vs time 
r for constrained and unconstrained simulations. The constrained simulation 
maintains K — 1 (except at bounces) while the unconstrained one evolves 
to K < 1. With K < 1, the system point is moving too slowly in MSS 
to bounce off the potential. Thus the unconstrained simulation yields the 
spurious result that the model has a last bounce and is AVTD. The simu- 
lation does not become unstable. There is no code crash — just a physically 



K = vl + v 2 _ = 1 



(5) 




u n > 2 
1< u n < 2 



(6) 
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Figure 1: Constraint enforcement in collapsing Bianchi IX cosmologies. Left: 
Comparison of constrained and unconstrained evolutions showing K vs t. Right: 
Sequence of ti-values for constraint enforcement every N time steps. 



wrong result. (From some initial data, the unconstrained evolution leads to 
K > 1. This does become unstable and crashes.) Figure 1 also compares 
the evolution from the same initial data while solving the constraint at ev- 
ery N time steps for N — 1, 2, 5, 10, 50. The evolution is monitored by 
the sequence of u-values obtained for the Kasner segments. Two evolutions 
become qualitatively distinct when the integer parts of u differ for the same 
segment. The sequence of u- values obtained from (6) is sensitive to initial 
conditions when an era ends [u — *■ 1/("U — 1)]. Note that the evolutions are 
qualitatively identical for a large number of bounces independent of iV and 
that they begin to differ at about the same place. This epoch is essentially 
the point at which information about the initial value of u is lost (in the 
double precision that is used). Thus one cannot really choose one of the 
solutions as more correct than another. The lesson here appears to be that 
constraint solving — at least from time to time — is essential for the accurate 
simulation of collapsing homogeneous Mixmaster models. 
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2.2 £7(1)-Symmetric Cosmology 



As our next example, we consider the vacuum L^l) -symmetric cosmologies 
discussed extensively by Moncrief [12]. The metric for T 3 spatial topology is 
given by (for a specific choice of lapse and shift) 

-e 2A dr 2 + e A e ab (x, z)dCd£ b } + e 2ip (d^ + [3 a dx a dr) 2 (7) 



ds 2 



-2ip 



where a, b = 1, 2 and <p, A, x, z, and f3 a depend on spatial variables £i, £ 2) 
and time r. The explicit form of e ab is given in [12, 15] as is the discussion of a 
canonical transformation to replace the twists (3 a with a single twist potential 
uj. Einstein's equations are obtained from variation of the Hamiltonian [13] 



H 



T 3 



V 



1 4b 2 1 2 

-e v r - -p A 



2 2 1 

+ {(e A e ab ) , ab -(e A e ab ) , a A, b +e A [(e- 2z ) , u x ru - (e~ 2z ) :V x n 



+2e A e 



A„ab 



<P,a<P,b+2 e e 



A.p-4tp p ab 



e LJ, a UJ, b 



n 



T 3 



(8) 



where Ti = is the Hamiltonian constraint and p 9 , r, p z , and p x are 
cannonically conjugate to tp, uj, A, z, and x. 

Polarized [/(l)-symmetric models (uj = = r) have been examined both 
numerically [15] and analytically [16] where the case has been made that the 
singularity is AVTD. Generic vacuum models are, so far, beyond analytic 
study. Numerical simulations indicate that [/(l)-symmetric collapse exhibits 
LMD [13, 17]. This behavior can be understoon in terms of the Method of 
Consistent Potentials (MCP) [7] or by analogy with spatially homogeneous 
models [18]. Snapshots from a typical evolution are shown in [13]. As was 
mentioned in [13], the LMD behavior cannot be achieved without explicit 
enforcement of the Hamiltonian constraint. If the momentum constraint is 
not enforced, the Hamiltonian constraint may be solved algebraically (e.g. for 
Pa) ■ Figure 2 illustrates that the wave amplitude (p exhibits LMD oscillations 
at a representative spatial point in a constrained evolution while, from the 
same initial date at the same spatial point, the oscillations are absent in 
unconstrained evolution. Figure 2 also compares the analog of the Kasner 
measure K (defined so that K = (p/^Pk) 2 + (Pz/^Pa) 2 = 1 is the analog 
of the Kasner circle) vs r for the same simulations. Once again, there is no 
instability in either case. As in the spatially homogeneous example, failure to 
enforce the Hamiltonian constraint leads to qualitatively incorrect behavior — 
a spurious indication that the model is AVTD. 
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Figure 2: Comparison of ip vs r (left) and K vs r (right) for constrained and 
unconstrained evolutions in ?7(l)-symmetric collapse from the same initial data. 
Behavior at a representative spatial point is shown. 



2.3 Generic T 2 -Symmetric Collapse 

Our final example will be generic T 2 -symmetric collapse [8]. These space- 
times are the most general T 2 -symmetric vacuum spacetimes and reduce to 
the Gowdy model with T 3 spatial topology [19] if the twists are set to zero. 
As in [8], consider the metric 



ds 2 



-e^^dr 2 + e^ +r)/2 d6 2 + 



2 ' e- p - 



T [d5 



dr'Q)def 



+e p ~ T \ da + QdS + 



dr' (Q9) -Q dr'Q 



dB 



(9) 



where the wave amplitudes P and Q and the "background" A depend only 
on the spatial variable 9 and time r, O = /te ( A + 2P + 3r )/ 2 e ^/ 4 for k the twist 
constant, and e^/ 4 = 2ix\. (For details see [8].) Einstein's equations may 
be derived from the Hamiltonian density {which is not the Hamiltonian con- 
straint) 



H 



4tt> 



Til 



e- 2P * 2 Q + e"* 



P'- 



aK 2 n x e^ +2P+3 ^ 2 . 



(10) 



This is supplemented by the momentum constraint 7T P P, e +7tQQ,g +7r A A )6 ) : 



0. The Gowdy model is recovered if k — 0, ix\ 



The first order equation 
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Figure 3: Comparison of P vs r (left) and u> vs r (right) for constrained and 
unconstrained evolutions in generic T 2 -symmetric collapse from the same initial 
data. Behavior at a representative spatial point is shown. 

for A obtained from the variation of (10) is in fact the Hamiltonian constraint. 
Thus these variables and equations "automatically" enforce the Hamiltonian 
constraint. A typical evolution is shown in [8]. 

The set of equations obtained from the variation of (10) plus the mo- 
mentum constraint do not comprise the full set of Einstein equations ob- 
tained from (9) and expressed as G^ u = 0. One also obtains a redundant 
wave equation for A (see [11]). This wave equation may be used instead of 
the Hamiltonian constraint to evolve A to produce an unconstrained evolu- 
tion. Figure 3 shows respectively the wave amplitude P and f-like parameter 
w = iip/27i\ (see [8]) at the same representative spatial point for constrained 
and unconstrained evolution from the same initial data. It is clear that 
the LMD discussed in detail and understood quantitatively in [8] through 
"bounce laws" for w cannot be reproduced without explicit enforcement of 
the Hamiltonian constraint. Once again, there is no numerical instability 
associated with failure to enforce the constraint. 
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3 Conclusions 



Three examples have been presented where failure to enforce the Hamiltonian 
constraint has led, not to instability, but to qualitatively incorrect behav- 
ior. The Hamiltonian constraint consists of kinetic energy-like and potential 
energy-like terms. In collapsing cosmologies with LMD behavior, the ki- 
netic terms either dominate or intermittently dominate the dynamics. If the 
Hamiltonian constraint is not enforced, the relationship among the momenta 
in the kinetic term will become incorrect. However, these momenta comprise 
the coefficients a of the time r in potential terms of the form e aT . How these 
terms grow and decay determine whether the approach to the singularity is 
AVTD or LMD (see the discussion in [7]). With the wrong coefficient — or 
even the wrong sign for a — a term that should grow exponentially may decay 
and vice versa in a numerical simulation. 

There are collapsing cosmological spacetimes that do not require explicit 
enforcement of the Hamiltonian constraint. Examples already mentioned in- 
clude the polarized U(l) models and the Gowdy cosmologies. Thus it seems 
as if the "delicate" LMD behavior must be present for enforcement of the 
Hamiltonian constraint to be needed. Enforcement of the momentum con- 
straint appears to be less critical. However, momentum constraint violation 
can be shown to lead to incorrect spatial waveforms in collapsing cosmologies. 
Garfinkle has demonstrated LMD in generic collapse without explicit con- 
straint enforcement (although damping was introduced to suppress constraint 
violating modes) [9] . The variables he used are more naturally adapted to the 
constraint hypersurface but are not suitable for binary black hole evolution 

[20] - 

The lesson from the examples presented here is that failure to enforce the 
Hamiltonian constraint can yield qualitatively incorrect behavior in a numer- 
ical simulation even in the absence of instability. While this does not occur 
in all classes of spacetimes and for all choices of variables and formulations, 
failure to enforce the constraint may be dangerous when properties of the 
solution are unknown. 
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